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Abstract 

Electric  fields  and  potentials  of  an  equilibrated  assembly  of  ions  and  water  molecules  adjacent 
to  a  charged  metal  surface  are  calculated  as  a  function  of  perpendicular  distance  z  from  the 
surface  from  data  derived  from  molecular  dynamics  trajectories.  The  spatial  distributions  of 
atoms  or  molecules  along  direction  z  are  found  by  ensemble  averaging  of  trajectories  followed 
by  averaging  with  a  localized  function  with  a  well  defined  length  scale.  Two  methods  were 
used  calculate  z  dependent  charge  density  distributions.  In  the  first,  to  be  called  the  atom 
method,  the  trajectories  of  charged  atoms  are  averaged.  In  the  second,  called  the  molecule 
method,  a  Taylor  expansion  of  charged  atom  positions  relative  to  molecular  centers  is  per¬ 
formed  and  the  charge  density  separated  into  monopole,  dipole,  quadrupole,  octopole,... 
components.  These  distributions  are  used  to  calculate  the  electric  potential  and  in  one  ex¬ 
ample  to  study  the  progressive  loss  of  structure  due  to  water  as  the  length  parameter  is 
scanned  through  the  dimension  of  a  water  molecule.  This  latter  result  provides  a  link  be¬ 
tween  simulations  with  detailed  atomic  modeling  of  intermolecular  interactions  and  electric 
potentials  derived  from  Gouy-Chapman  theory.  Illustrative  examples  are  chosen  from  sim¬ 
ulations  of  aqueous  solutions  of  simple  alkali  halide  electrolytes  next  to  charged  and  un¬ 
charged  flat  metal  surfaces.  The  smallest  system  has  one  ion  and  157  water  molecules,  the 
largest  60  ions  and  1576  water  molecules. 
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Abstract 

Electric  fields  and  potentials  of  an  equilibrated  assembly  of  ions  and  water  molecules  adjacent 
to  a  charged  metal  surface  are  calculated  as  a  function  of  perpendicular  distance  z  from  the 
surface  from  data  derived  from  molecular  dynamics  trajectories.  The  spatial  distributions  of 
atoms  or  molecules  along  direction  z  are  found  by  ensemble  averaging  of  trajectories  followed 
by  averaging  with  a  localized  function  with  a  well  defined  length  scale.  Two  methods  were  used 
calculate  z  dependent  charge  density  distributions.  In  the  first,  to  be  called  the  atom  method, 
the  trajectories  of  charged  atoms  arc  averaged.  In  the  second,  called  the  molecule  method,  a 
Taylor  expansion  of  charged  atom  positions  relative  to  molecular  centers  is  performed  and  the 
charge  density  separated  into  monopole,  dipole,  quadrupole,  octopole,...  components.  These 
distributions  are  used  to  calculate  the  electric  potential  and  in  one  example  to  study  the  pro¬ 
gressive  loss  of  structure  due  to  water  as  the  length  parameter  is  scanned  through  the  dimension 
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of  a  water  molecule.  This  latter  result  provides  a  link  between  simulations  with  detailed  atomic 
modeling  of  intermolecular  interactions  and  electric  potentials  derived  from  Gouy-Chapman 
theory.  Illustrative  examples  are  chosen  from  simulations  of  aqueous  solutions  of  simple  alkali 
halide  electrolytes  next  to  charged  and  uncharged  flat  metal  surfaces.  The  smallest  system  has 
one  ion  and  157  water  molecules,  the  largest  60  ions  and  1576  water  molecules. 


I.  INTRODUCTION 


In  this  paper  we  describe  the  calculation  of  electrostatic  fields  arising  from  ions  and  polar  mol¬ 
ecules  in  equilibrium  distributions  in  front  of  charged  metal  electrodes.  The  focus  is  on  fields 
within  ten  water  molecule  diameters  of  the  electrode-aqueous  electrolyte  interface.  Molecular 
d5Tiamics  simulations  are  used  to  calculate  the  trajectories  of  ions  and  polar  molecules,  which 
are  time  and  space  averaged  in  a  manner  to  be  described  to  get  time  independent  distributions 
depending  only  on  the  coordinate  z  perpendicular  to  the  surface.  These  averaged  distributions 
are  then  used  as  the  source  terms  in  Maxwell's  equations  to  calculate  fields  and  potentials  that 
are  in  a  form  that  can  be  compared  with  potentials  derived  intuitively  from  electrochemical 
knowledge  or  the  results  of  non  simulation  methods  like  that  of  Gouy-Chapman-Stem^'^  theory 
or  Henderson  and  coworkers  more  sophisticated  correlation  function  techmque^.  These  latter 
methods  solve  Poisson's  equation  with  the  Boltzmann  ansatz  and  calculate  from  the  charge  dis¬ 
tribution  the  electric  field  and  potential  across  the  double  layer  region.  In  our  work  examples 
are  chosen  from  simulations  of  aqueous  solutions  of  simple  electrolytes  next  to  charged  metal 
surfaces  that  represent  frequently  encountered  situations  in  adsorption  from  electric  double  lay¬ 
ers.  All  sums  of  electrostatic  interactions  are  evaluated  without  spatial  cut-offs  using  the  fast 
multipole  method  of  Greengard  and  coworkers^'^ 

The  cartoon  in  Figure  1  depicts  the  traditional  view  of  the  electric  double  layer  in  the 
thermodynamically  stable  region  of  electrode  potential  for  a  metal  with  a  flat  surface.  The 
electric  potential  due  to  ions  is  schematically  shown  at  top.  We  will  show  that  near  the  surface 
this  way  of  depicting  the  distance  dependence  of  the  potential  is  wrong  because  it  neglects  the 
contribution  from  oriented  water  molecules.  In  Figure  1  the  labels  IHP  and  OHP  denote  the 
inner  and  outer  Helmholtz  planes  respectively.  In  the  cartoon  the  compact  part  of  the  double 
layer  contains  an  anion  in  physical  contact  with  a  flat  surface  and  a  cation  two  water  molecules 
distant  from  the  electrode  as  in  tlie  the  model  of  Bockris,  Devanatlian  and  Muller  .  The  diffuse 


part  is  symbolized  by  the  labelled  arrows  which  cover  the  region  from  the  OHP  to  the  bulk 
electrolyte  region. 

The  calculation  of  the  classical  electric  field  of  a  set  of  charges  is  a  well  studied  problem.  Given 
a  set  of  charges  the  electric  field  and  potential  can  be  found  with  arbitary  high  accuracy  by  a 
number  of  methods.  In  molecular  dynamics  and  Monte  Carlo  simulations  Ewald's  method  has 
been  frequently  used.  However  in  systems  large  numbers  (N  >  200)  of  charged  or  polar  species, 
Ewald  methods  are  slow  (time  t  compared  to  the  fast  multipole  method  (t  «  AT)  and 

particle-mesh  methods^  ^  (t «  log  N).  The  molecular  dynamics  simulations  described  in  this 
paper  would  not  be  possible  without  the  use  of  the  fast  multipole  method  to  evaluate  the  long 
range  electrostatic  fields.  Using  the  fast  multipole  method  we  have  performed  many  simulations 
on  systems  containing  600  charged  particles  and  a  few  simulations  on  larger  systems  containing 
5000  charged  particles  and  their  electrostatic  images. 

Given  that  we  are  able  to  accurately  evaluate  long  range  electrostatic  interactions  in  a  molecular 
dynamics  simulation,  the  next  step  is  the  calculation  of  local  electric  fields  from  time  averaged 
charge  distributions.  This  is  the  main  task  of  the  present  paper.  We  concentrate  on  two  methods 
to  calculate  time  and  space  averages  of  charge  source  functions,  and  the  corresponding  electric 
fields  and  potentials.  The  first  method,  to  be  called  the  atom  charge  method,  uses  the  positions 
of  the  atoms  and  the  charge  or  partial  charge  of  the  atoms  to  calculate  a  time  independent  electric 
charge  density.  After  making  averages  parallel  to  the  surface  this  technique  yields  results 
equivalent  to  that  described  by  Wilson,  Pohorille  and  Pratt'^’  based  on  the  derivation  given 
in  Landau  and  Lifshitz^®.  The  second  method  (called  the  molecule  method)  considers  the  system 
in  a  fundamentally  different  way.  It  views  the  world  as  a  collection  of  charged  or  polar  mole¬ 
cules.  It  uses  a  Taylor  series  expansion  of  charged  atom  coordinates  relative  to  an  origin  on  the 
molecule  to  express  the  electric  source  in  terms  of  a  sequence  of  electric  multipole  (monopole, 
dipole,  quadrupole,  octopole,...)  polarization  source  densities.  This  method  follows  the  de- 
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scription  given  in  standard  texts  such  as  Jackson^^,  which  is  based  in  part  on  the  lucid  de¬ 
scription  in  the  paper  by  Russakoff^^.  By  this  second  method  the  contribution  to  the  electric 
field  and  potential  of  monopole,  dipole,  quadrupole,  octopole,...  moments  of  the  molecules  can 
be  calculated  separately  in  a  systematic  way,  and  could  in  principle  be  used  to  input  exper¬ 
imentally  measured  values  of  electric  multipole  moments.  Each  of  the  two  methods  provides 
its  own  particular  insight  in  the  computation  of  the  fields. 

To  illustrate  the  calculation  of  the  fields  and  potentials  we  choose  examples  from  on  going  work 
in  simulations  of  aqueous  electrochemical  interfaces  that  typify  general  classes  of  behaviour. 
All  the  systems  were  composed  of  water  and  monovalent  ions  adjacent  to  a  metal  surface.  The 
examples  chosen  from  small  scale  simulations  are:  a  non  contact  adsorbed  Li  ion,  a  contact 
adsorbed  I  ion,  a  neutral  solution  of  Li  and  I  (amon  contact  adsorbed),  a  neutral  solution  of 
Li^  and  F  ,  and  Li''  and  two  iodides  21  (cation  attracted  to  contact  adsorbed  anions).  As  an 
example  of  a  larger  calculation  that  shows  that  some  of  the  main  features  in  the  small  simulations 
are  also  found  in  calculations  ten  times  larger  is  given  at  the  end  of  the  paper.  In  this  example 
we  discuss  a  IM  NaCl  salt  solution  in  front  of  a  charged  electrode.  In  this  simulation  the  region 
of  electrolyte  that  screens  the  charge  on  the  electrode  and  the  bulk  solution  are  evident 

In  our  smallest  simulations  the  cell  size  (L  =  1.862  nm)  is  too  small  for  us  to  be  confident  that 
we  can  do  more  than  describe  the  average  electric  field  between  the  metal  and  the  OHP.  The 
decay  of  the  average  electric  field  to  zero  at  the  non  metal  boundary  is  the  result  of  charge 
neutrality  in  the  immersed  electrode  model.  Consequently  in  the  smaller  simulations  we  focus 
primarily  on  the  interfacial  fields  within  the  first  three  or  four  water  molecules,  corresponding 
to  a  distance  up  to  about  one  mn  from  the  electrode.  In  the  large  simulation  with  a  cell  (L  — 
3.74  nm)  containing  5000  charges^^,  the  distribution  of  ions  from  the  zone  where  the  screening 
of  the  charged  metal  occurs,  all  the  way  to  the  'bulk'  electrolyte  is  observed.  In  this  case  we 
can  more  confidently  interpret  the  decay  of  field  and  potential  found  in  the  calculations  and 
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identify  what  corresponds  to  the  diffuse  region  of  the  electric  double  layer.  Traditionally  the 
structure  of  the  compact  region  where  ions  contact  adsorb  is  thought  of  as  being  static  and  re¬ 
sembling  a  capacitor  of  atomic  separations.  Our  model  permits  the  ions  to  move  in  and  out 
freely  from  the  diffuse  layer  and  to  contact  adsorb  and  desorb.  We  note  that  Henderson  and 
coworkers^*^  in  their  work  solving  the  Omstein-Zemike  equation  for  the  electrode-electrolyte 
interface  have  shown  (as  does  the  work  presented  here)  that  there  is  no  sharp  division  between 
the  diffuse  and  inner  layers,  in  contrast  to  the  ideas  conveyed  by  the  pictures  in  many  textbooks. 

We  summary  the  goals  of  the  present  woric  are  as  follows.  First,  the  calculation  of  averaged 
electric  field  and  potential  for  systems  of  electrochemical  interest  Second,  to  contrast  the  atom 
and  molecule  methods  of  calculating  charge  densities  and  potentials.  Third,  to  qualitatively 
analyse  deviations  in  the  electric  potential  from  the  curve  for  ions  in  a  dielectric  fluid  interms 
of  distribution  of  solvent  and  ions  near  the  metal  surface. 

We  close  this  section  with  a  brief  mention  of  important  parts  of  the  physics  that  are  still  missing. 
First  there  is  no  feature  that  permits  explicit  consideration  of  the  metal  surface  topography  and 
conduction  electron  densities  at  the  surface.  For  example,  one  can  go  much  further  and  use 
jellium  models  or  jellium  with  embedded  ions  but  this  is  beyound  the  scope  of  the  current  study. 
Halley  and  coworkers^^'^^  have  developed  theories  based  on  jellium  electrode  models  of  the 
charged  metal  surface  and  adjacent  solution  that  accounts  for  aspects  of  measured  capacitance 
vs.  potential  curves.  Cluster  calculations  have  been  used  to  obtain  metal-water  and  metal-ion 
potentials  for  use  in  molecular  dynamics  simulations  that  introduce  aspects  of  metal  surface 
corrugation^"^'^^  30-34  35-37  no  effects  due  to  electronic  polarizability 

or  geometric  distortions  of  molecules  in  the  high  electric  fields  of  the  double  layer. 
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11.  MODELS  AND  METHODS. 

A.  The  Immersed  Electrode  Model. 

When  the  electrode  potential  is  altered  and  charge  flows  onto  the  electrode,  the  composition  of 
the  electrolyte  next  to  the  electrode  adjusts  to  screen  the  new  charge  on  the  electrode.  For  dilute 
solutions  (<  0.1  M)  a  rough  estimate  of  the  screening  layer  thickness  is  ic‘,  where  k  is  the 
Debye-Huckel  screening  constant.  According  to  the  Gouy-Chapman  theory  the  concentration 
of  monovalent  ions  has  fallen  to  e'^  of  its  'surface'  value  at  a  distance  d  from  the  electrode  given 
by 


\  V  / 

Here  £  is  the  macroscopic  dielectric  constant  of  the  solvent,  and  n^  is  the  concentration  of  the 
ions  in  the  bulk.  Typical  values  of  d  are:  3.1  nm  for  0,01  M,  and  0.96  nm  for  0.1  M.  At  higher 
salt  concentrations  outside  the  range  of  the  derivation  of  this  formula  we  get  0.55  nm  for  0.3 
M,  and  0.31  nm  for  1  M  salt  solutions.  The  formula  provides  some  rough  measure  perhaps  even 
within  a  factor  2  (tlie  effective  dielectric  constant  of  water  near  the  electrode  may  be  ten  times 
smaller  than  in  bulk^^)  of  the  double  layer  thickness  even  though  the  basis  of  the  derivation  is 
not  valid.  These  brief  considerations  suggest  we  can  use  molecular  dynamics  to  simulate  an 
immersed  electrode  when  the  salt  concentration  is  about  0.3  M  or  larger  by  including  the  adja¬ 
cent  electrolyte  region  out  to  a  tliickness  of  about  ten  water  molecules.  In  this  paper  we  do  this 
first  for  a  small  cell  about  five  water  molecules  tliick,  expecting  only  to  be  able  to  model  the 
so-called  inner  part  of  the  double  layer  where  contact  adsorption  occurs.  Later  we  study  a  cell 
about  10  waters  thick  and  are  able  to  identify  a  region  in  space  that  has  bulk-like  properties, 
lending  support  to  the  immersed  electrode  model  as  we  use  it.  Our  immersed  electrode  model 


7 


consists  of  a  layer  of  electrolyte  between  two  walls.  The  wall  on  the  left  carries  no  charge  it 
is  simply  a  restraining  wall  and  ideally  would  allow  a  continuous  transition  to  the  bulk 
electrolyte  region.  The  complete  system  of  electrolyte  and  electrode  (always  on  the  right  hand 
side  in  all  the  Figures  of  this  paper)  is  neutral.  The  approach  is  useful  because  it  reduces  the 
number  of  water  molecules  in  the  calculation,  and  because  there  is  only  one  metal  surface,  there 
is  only  one  electrostatic  image  plane.  Finally  we  point  out  that  in  our  system  the  electrical 
properties  of  the  system  are  due  to  three  layers.  The  vacuum  half  space  occupying:  -°o  <  z  <  - 
VzL;  the  xy  periodically  replicated  simulation  cells:  ML<  z  <  16L;  and  the  metal  half  space:  ViL 
<  z  <  oo.  More  details  and  discussion  of  the  immersed  electrode  model  are  found  elsewhere  . 

B.  Models  for  water  and  ions. 

In  the  smaller  scale  simulations  reported  here  we  used  the  parameters  of  the  Stillinger  ST2 
water  model  and  the  interaction  parameters  with  alkali  metal  ions  and  halide  ions  developed  by 
Heinzinger  and  cowoikers  The  ST2  water  molecule  model  consists  of  a  central  oxygen  atom 
surrounded  by  two  hydrogen  atoms  and  two  massless  point  charges  (PC)  in  a  rigid  tetrahedral- 
like  arrangement.  The  charges  on  the  particles  are  q  =0.23570lel,  and  q  =0.  The 

alkali  metal  and  halide  ions  are  non-polarizable  Lennard-Jones  atoms  with  point  mass  and 
charge.  The  atom-atom  interaction  parameters  are  taken  from  Heinzinger's  review  In  the 
larger  simulations  with  a  3.74  nm  wide  cell  containing  IM  NaCl  solution  tlie  SPCE  water  model 
was  used  primarily  because  this  water  model  has  fewer  charges  (tliree  compared  to  four  for  ST2) 
and  requires  less  time  to  compute  the  electrostatics.  The  Lennard-Jones  potential  between  the 
atoms  of  each  molecule  was  smoothly  cut-off  at  distances  of  R  =  0.68  nm. 

C.  Molecule-Wall  and  Ion-Wall  Potentials. 

There  are  two  walls.  The  simplest  is  the  uncharged  restraining  non  metallic  wall  on  the  left  hand 
side  in  Figures  2  -  15.  It  restrains  the  fluid  by  a  9-3  potential.  The  metal  wall,  on  the  right  hand 
side,  is  represented  by  two  superimposed  potentials.  The  first  is  a  9  -  3  potential  to  represent 
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both  the  Pauli  repulsion  and  dispersive  attractive  interactions.  The  second  potential  is  an 

electrostatic  image  potential  that  describes  the  interaction  between  a  charge  in  the  electrol5de  and 

the  conduction  electrons  of  the  metal.  In  the  calculations  described  here  the  image  plane  and 

origin  plane  of  the  9-3  potential  were  coincident  This  is  equivalent  to  choosing  the  image  plane 

and  the  nuclear  plane  of  the  metal  surface  to  be  the  same.  This  is  acceptable  in  our  scheme 

because  the  Leonard- Jones  core  parameters  o  are  all  large  and  the  'thickness  of  the  repulsive 

wall  is  also  large  (ca.  0.247  nm).  The  atom-surface  interaction  parameters  describing  interaction 

42 

with  nonconduction  electrons  were  chosen  to  be  the  same  as  those  of  Lee  at  al  , 
A=17.447xl0'^  kJCnm^Anole  and  B=76.144xl0'^  kJ(nm)Vmole  for  O,  I  ion  and  Li  ion.  The  A 
and  B  parameters  for  H  and  PC  were  set  to  zero  for  water  molecules.  Note  that  in  this  model 
the  adsorption  well  depths  of  the  Lennard-Jones  wall  potentials  is  a  few  kcalAnole,  similar  to 
kT  at  room  temperature,  and  therefore  pretty  whimpy. 

D.  Advantages  of  the  Fast  Multipole  Method, 

Simulations  of  aqueous  electrolyte  solutions  using  molecular  dynamics  and  Monte  Carlo  meth¬ 
ods  require  the  evaluation  of  the  superimposed  Coulomb  fields  of  thousands  particles  (a  5  nm 
cube  of  water  contains  about  12,500  atoms).  In  this  paper  we  use  the  fast  multipole  method 
(fmm)  of  Greengard  and  Rokhlin^'^  to  evaluate  the  electric  fields  acting  on  the  charged  particles 
during  the  computer  simulation.  A  number  of  methods  are  in  use  to  calculate  or  approximate 
long  range  Coulomb  fields.  Friedmann  and  Honig"^^  have  surveyed  some  of  empirical  dielectric 

.44 

recipes  in  use  in  biological  simulations.  For  example  recipes  like  the  Hingerty  function  and 
other  distance  dependent  dielectric  'constants’  have  been  used  in  simulations  of  proteins. 
Direct  summation  with  a  cut  off  after  about  1.00  nm  is  common  in  many  commercially  available 
codes  like  Charmm^^  and  Amber^^’  There  are  a  number  of  plane  wise  summation  methods 
from  crystal  physics^^.  For  homogeneous  systems  the  reaction  field  metliod  of  Barker  is 
simple  and  easy  to  use.  This  method  is  not  easily  applied  to  liquid-solid  interfaces  though  at- 
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tempts  to  extend  it  have  been  described  The  Ewald  method  has  been  extensively 

used  when  rigor  and  accuracy  are  needed.  Ewald's  summation  crudely  applied  is  proportional 
to  and  at  best  where  N  is  the  number  of  charges.  The  fast  multipole  method  (fmm) 
developed  by  Greengaixi  and  Rokhlin  is  an  order  N  algorithm,  and  consequently  is  the  only 
viable  method  for  very  large  simulations.  The  fmm  technique  is  attractive  because  of  the  ease 
of  implementation  of  a  variety  of  boundary  conditions  such  as  periodic,  Dirichlet,  Neumarm  and 
mixed  boundaries  and  because  an  adaptive  version  of  the  algorithm  is  available  in  which  re¬ 
gions  of  low  or  no  charge  density  are  not  subdivided  when  the  charge  count  falls  below  a 
specified  integer. 
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III.  ELECTROSTATICS  OF  ELECTROCHEMICAL  CELLS 

There  are  many  ways  that  electric  fields  and  potentials  can  be  calculated.  For  example,  at  every 
time  step  we  use  the  fast  multipole  method  to  compute  the  electric  forces  acting  on  every 
charged  particle  in  the  system.  We  could  take  the  electric  field  evaluated  at  a  specified  particle, 
say  an  O  atom,  at  a  given  instant  of  time  and  then  either  allow  the  molecule  to  move  and  average 
the  field  at  regular  intervals  or  confine  the  atom  to  a  spatial  well  that  would  then  be  transported 
so  as  to  sample  all  space.  This  would  clearly  be  useful  for  interpreting  experiments  like  NMR, 
but  is  not  in  the  spirit  of  Gouy-Chapman  theory  which  attempts  to  solve  the  Poisson  equation 
using  a  Maxwell-Boltzraann  ansatz  for  local  charge  fluctuations.  It  is  not  pursued  further  here. 
Instead  in  this  section  we  describe  two  methods  for  getting  source  distributions  to  be  used  to 
calculate  averaged  electric  fields  and  potentials  for  comparison  with  other  theoretical  methods 
which  are  in  the  spirit  of  Gouy-Chapman  approach.  The  first  method  assumes  that  the  atoms 
are  fundamental  objects  and  uses  the  charge  or  partial  charge  on  each  atom  and  the  three  space 
coordinates  of  each  atom  recorded  at  regular  time  separations  (usually  0.5  or  1  ps)  as  input  into 
the  generation  of  a  charged  source  distribution.  We  call  this  the  atom  approach.  The  only  source 
term  is  the  charge  density  function  in  vacuum.  Without  any  local  spatial  averaging  this  point 
of  view  is  similar  to  that  followed  by  Wilson,  Pohorille,  and  Pratt^^’  The  second  way,  to 
be  called  the  molecule  method,  assumes  that  the  system  is  a  collection  of  molecules  that  have 
internal  electrical  structure.  These  are  inherently  more  complex  objects  than  atoms.  To  specify 
the  electrical  properties  we  need  space  coordinates,  orientations  and  electrostatic  multipole  mo¬ 
ments  for  each  molecule  in  the  system.  The  approach  we  follow  can  be  found  in  the  clearly 
written  article  by  Russakoff'*^,  or  the  treatise  by  de  Groot  and  Suttorp"’*’,  botli  sets  of  authors 
derive  a  set  of  macroscopic  equations  from  Maxwell's  equations  in  vacuum  for  point  charged 
particles.  Their  approach  is  summarized  in  a  number  of  standard  texts  .  We  note  that  though 
the  dynamics  are  uniquely  defined  by  the  models  we  use,  the  time  independent  average  fields 
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and  potentials  that  are  calculated  are  not  unique  because  we  can  broaden  the  point  charge  dis¬ 
tributions  inside  the  ions  and  molecules  in  a  spherically  symetric  way  and  so  long  as  the  broaden 
distributions  stay  well  inside  the  Lennard-Jones  spheres  (so  that  distributions  on  different  mol¬ 
ecules  do  not  overlap)  the  the  dynamics  are  unchanged.  This  means  we  can  smooth  some  of  the 
averages  within  limits  set  by  geometry  of  the  models  and  the  energetics  of  collisions.  This  point 
was  clearly  recognized  by  Wilson,  Pohorille  and  Pratt^"^’  in  their  interesting  discussion  of  the 
properties  of  the  vacuum-water  surface. 

The  atom  method  uses  the  distribution  of  point  charges  in  vacuum.  Consider  the  set  of  point 
charges  without  regard  for  whether  they  originate  from  neutral  water  molecules  or  charged  ions. 
In  this  case  the  source  term  for  Maxwell  fields  in  vacuum  is  the  microscopic  charge  density 


Here 


and 


Patomsir^t)  =  ^  qj5{r  -  r/o'),  [3.3] 

y=i 

where  is  the  number  of  atoms  in  the  simulation  cell  and  ry  is  the  position  of  the  j-th  atom. 
The  surface  charge  density  on  the  metal  p^Ux,y,t)  is  time  dependent  because  it  depends  on  the 
position  of  tlie  atomic  charges.  All  die  charge  in  the  system  is  described  by  p;. 
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We  can  subject  this  charge  density  to  local  space  averaging  using  a  test  function.  Let  fir)  be  a 
real  positive  function  localized  around  r  =  0.  We  define  the  local  spatial  average  of  F{r,  t)  by 


<F(r,0>  =  jrfr'y(r')/'(r-r',0. 


The  time  average 


_  lim 
F(r)  = 


r 


r')>  df 


can  be  replaced  by  an  average  over  configurations  at  different  times 


[3.4] 


[3.5] 


N, 


F(r)= 


1 


cor^igs 


N. 


corfigs 


^  <F(r,  ti)>. 


i=  1 


[3.6] 


The  systems  we  consider  have  translational  invariance  in  the  xy  plane.  The  function  F(r)  is  a 
function  of  z  only.  In  light  of  this  we  consider  only  test  functions  which  are  functions  of  z. 
We  use  two  test  functions,  a  bin-like  test  function 


r  1/g  ,  izi  <  '4? 

\  0  ,  \z\  >  Vig 


[3.7] 
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and  a  localized  one  dimensional  Gaussian  function 


A^)  =  iKgY'^e-^<  [3.8] 

The  metal  surface  charge  density  p^UJc.y,/)  is  replaced  by  the  averaged  image  charge  density 
which  is  a  constant.  The  metal  image  charge  is  denoted  by  After  we  have  averaged  over 
the  many  spatial  configurations  from  the  molecular  dynamics  calculation  we  obtain  a  z  de¬ 
pendent  charge  density  profile  given  by 


PK^)  ~  Prnetal^^^  VlL)  +  Patomsi^^ 


where  as  in  Eq.(3.6) 

^configs 

Paiomsi^^  “  ^  \  ^Patomsi^^ 

^^corftgs 

After  substituting  explicitly  for  the  test  function  we  get 

N 

atoms 

^(7y/(z-Zy(r,)) 

>=1 


K 


Patomsi^) 


1 


A 


configs 


configs 

I 

1  =  I 


[3.10] 


[3.11] 


The  electric  field  is  given  directly  by  integrating  the  vacuum  Maxwell's  equation.  Since  we  have 
already  averaged  over  the  xy  plane  the  electric  field  nonnal  to  the  surface  is  given  by 
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E,/z)  =  ^\  dz'-piiz') 

^  oo 


[3.12] 


The  electric  potential  is  given  by  a  second  integration 


[3.13] 


The  second  approach  views  the  system  as  a  collection  of  molecules  (labels  n,  positions  r„) 
composed  of  atoms  (label  b,  position  T;,*,  charge  qnb)-  The  charge  density  of  the  system  is. 


p//(r,0  =  Pmetal(X,y,t)^Z  -  ViL)  +  p„o/(r,0 


[3.14] 


where  Pm^ui  is  the  same  charge  density  on  the  metal  surface  as  in  the  atomic  method,  and 


P,„./(r4)  =  ^P„(r,0. 

n  =  ] 

Here  N,„oi  is  Uic  total  number  of  molecules  in  tlie  system. 


p„(r,/)  =  ^f7nfeS(r  -  r„  -  r,,^), 

b=\ 


[3.15] 


[3.16] 
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is  the  charge  density  of  the  n  th  molecule.  The  number  of  atoms  in  the  n  th  molecule  is  N„a, 
and  is  the  position  of  the  charge  q„b  measured  from  the  center  r„  of  the  n  th  molecule. 

Next  we  perform  explicit  local  spatial  averaging  with  a  test  function  and  take  the  ensemble  av¬ 
erage.  These  steps  are  formally  the  same  as  described  for  the  atom  method.  For  each  molecule 
label  n  we  take  the  average  with  respect  to  the  test  function /(z),  then  make  a  Taylor  expansion 
of  atomic  coordinates  relative  to  the  molecular  centers  r„.  This  yields  the  following  ex¬ 
pression  for  the  averaged  charge  density  by  the  molecule  method 


P«.o/(z)  =  P(2) 


-±p^(z)+-^ 

dz  dz^ 


Qjz)-^0,Jz)  + 
dz 


[3.17] 


where 


p(z)  = 


1 


N, 


configs 


1 

;  =  1  nt=\ 


[3.18] 


F,(z)  = 


1 

^configs 


^conftgs 


s 

\ 


^mol 

m  =  1  b= \ 


[3.19] 


K 


configs  ^mol 


e,,(z)=— ^ -  y  '/zY  f{2-zJti))yqnAhifiKbiti) 

/'nnfiov 


^configs  'TT 

^  °  f  =  1  in  -  \ 


b=  1 


[3.20] 
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^corfigs  ^mol 

=  J,  - -  y(z  —  Z^{ti)  )  ^(lmb^mbOi)^mbi^i)^mbih)  [3.21] 

‘^corfiss  ^ 

Here  the  charge  density  in  the  second  method  has  been  resolved  into  contributions  from 
monopoles  (ions  only  in  this  paper),  dipoles,  quadrupoles,  octopoles,  and  higher  order  terms. 
In  principle,  experimental  moments  can  be  substituted  for  dipoles,  quadrupoles,..  where  these 
are  known. 

In  the  following  sections  we  report  simulations  designed  to  typify  general  electrochemical  sys¬ 
tems.  The  distributions  displayed  in  most  of  the  Figures  were  calculated  using  a  one  dimensional 
binning  function  with  width  0.004655  nm.  In  one  case  to  be  described  in  next  section  we  discuss 
the  effect  of  changing  the  size  of  the  region  averaged  from  smaller  (.03  nm)  to  larger  than  the 
dimension  of  a  water  molecule. 
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IV.  ALKALI  METAL  ION 

This  section  describes  molecular  dynamics  calculations  with  one  Li"^  cation  and  157  ST2  water 
molecules  against  a  metal  surface.  The  simulation  cell  is  periodically  replicated  in  the  xy  plane 
parallel  to  the  metal  surface.  The  simulation  was  run  for  2500  ps,  the  first  100  ps  of  which  were 
used  to  equilibrate  the  system.  The  charge  on  the  electrode  is  the  image  charge  -lei.  At  any 
instant  this  charge  is  not  uniformily  distributed  across  the  electrode  but  localized  on  the  surface 
in  such  a  way  as  to  produce  the  same  electric  field  and  potential  as  the  electrostatic  image  of 
the  lithium  ion  and  all  the  water  molecules.  The  field  acting  on  the  lithium  ion  comes  from  all 
the  water  molecules  in  the  cell,  all  water  and  ions  in  xy  periodically  replicated  cells,  and  all  of 
the  electrostatic  images  of  the  contents  of  all  cells  in  the  image  plane  of  the  metal.  We 
emphasise  again  that  in  this  calculation  as  in  all  the  others  descibed  in  this  paper  no  electrostatic 
interaction  is  truncated.  Figures  2  and  3  show  the  probability  density  profiles  averaged  over  the 
xy  direction  and  electrical  potential  calculated  by  the  two  ways  described  in  section  III.  Note 
that  in  Figure  2  the  probability  of  finding  the  ion  and  water  have  scales  differing  by  fifty.  Also 
plotted  in  Figure  2  is  the  total  charge  density  by  the  atomic  method.  The  bin  size  is  L/400  = 
0.004655  ran  where  L  is  the  edge  length  of  the  simulation  cell. 

In  Figure  2  we  see  that  the  Li'^  ion  mass  center  mapped  out  a  diffuse-like  region  between  -0.6 
and  0.4  ran.  The  lithium  Li"^  ion  has  the  smallest  ionic  radius  of  all  the  monovalent  cations. 
Consequently  it's  hydration  shell  is  very  strongly  bound  making  it  more  difficult  for  this  ion  to 
contact  adsorb  on  the  electrode.  The  asymmetry  of  the  distribution  may  well  be  affected  by  the 
small  width  of  the  cell  and  consequently  no  significance  is  attached  to  its  shape  other  than  it  is 
diffuse  in  nature.  There  may  possibly  be  an  hydrophobic  plating  of  the  ion  to  tlie  left  hand 
interface.  On  the  metal  side  the  ion  rarely  approached  closer  than  0.3  ran  to  the  repulsive  part 
of  the  wall  potential  shown  by  dashed  line  at  0.68  ran  in  Figure  2.  This  is  approximately  the 
average  of  the  separate  a  parameters  for  Li^  and  0_ST2  suggesting  that  lithium  ion  remains  fully 
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hydrated  at  the  metal  surface,  and  that  though  surface  water  restricts  the  approach  of  the  ion  to 
the  metal  it  does  not  exclude  it.  This  behaviour  of  lithium  ion  is  more  like  that  postulated  in 
the  Grahame’s  picture  of  the  double  layer^^,  than  say  the  Bockris  model^^.  To  confirm  such 
behaviour  would  require  extending  the  simulation  over  many  nanoseconds,  or  performing  um¬ 
brella  sampling^^.  This  particular  aspect  of  the  problem  is  not  pursued  further  here. 

The  water  profile  shows  some  new  structure  not  seen  in  water  without  ions  .  Most  inter¬ 
esting  is  the  leading  peak  (closest  to  the  metal  surface)  at  ca.  0.68  run  due  to  a  few  localized 
water  molecules.  The  orientation  of  these  localized  water  molecules  can  be  ascertained  from  the 
H_ST2  and  PC_ST2  probability  distributions  (not  shown  here).  The  protons  on  these  molecules 
give  rise  the  distinct  peak  at  ca.  0.75  nm  in  the  H_ST2  distribution.  In  the  PC_ST2  distribution 
(not  shown)  there  is  peak  at  ca.  0.625  nm  that  is  enhanced  more  than  the  second  peak  in  the 
H_ST2  profile  which  is  at  ca.  0.6  nm.  The  first  water  peak  at  ca.  0.68  nm  lies  between  the  first 
H  and  PC  peaks  measured  from  the  metal  surface.  It  appears  therefore  that  some  of  the  localized 
water  molecules  have  one  proton  pointing  at  the  electrode. 

The  atomic  charge  density  profile  in  Figure  2  is  dominated  by  neutral  water.  The  large  oscil¬ 
lation  centered  near  0.70  nm  is  due  to  partially  oriented  water.  The  main  positive  peak  at  ca. 
0.75  mn  is  due  to  protons  attracted  to  the  metal  by  their  images  and  the  image  field  of  the  lithium 
ion.  There  are  smaller  oscillations  in  charge  profile  further  from  the  surface  that  are  more  clearly 
seen  if  the  distribution  is  smoothed.  These  tend  to  follow  the  oscillations  in  the  water  density 
probability  in  Figure  2.  The  actual  lithium  charge  density  is  buried  under  the  contribution  from 
the  water  when  the  bin  size  is  0.004655  nm. 

Figure  3  shows  the  electric  potential  across  tlie  cell  calculated  by  the  atomic  charge  method,  and 
by  the  molecular  method  with  systematic  inclusion  of  higher  electrostatic  multipoles.  For  clarity 
the  dashed  vertical  line  denoting  the  point  z  =  0.628  mu  where  tlie  wall  potential  goes  through 
zero  has  been  omitted  except  for  tlie  tick  mark  on  the  x  axis.  Since  all  the  curves  go  off  scale 
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near  z  =  0.682  nm  all  the  data  for  z  >  0.4  nm  have  been  plotted  again  on  a  larger  scale  and  are 
shown  in  the  right  hand  panel.  Note  too  that  in  this  right  hand  panel  the  monopole  curve  m 
(broken  curve)  is  scaled  by  0. 1  to  distinguish  it  from  the  other  curves.  The  notation  in  the  fig¬ 
ures  is  m  =  monopole,  d  =  dipole,  q  =  quadrupole.  When  m  and  d  contributions  are  combined 
the  curve  is  labelled  m+d.  When  m,  d  and  q  contributions  are  combined  the  curve  is  labelled 
m+d+q.  Only  when  three  multipoles  are  included  (m+d+q  curve)  do  the  potentials  calculated 
by  the  two  methods  agree  reasonably.  Due  to  the  small  size  of  the  bin  it  is  necessary  to  include 
more  higher  multipoles  before  the  atom  and  molecular  methods  agree  quantitatively.  The 
monopole  contribution  m  (dashed  line  in  Figure  3)  due  to  lithium  ion  alone  drops  monotonically 
as  expected  for  the  potential  inside  a  capacitor  where  the  charge  on  the  left  is  in  a  diffuse  layer 
spatially  separate  from  the  charge  on  the  right  plate  (metal)  at  z  =  0.931  nm.  Adding  the  dipole 
completely  changes  the  potential  (see  m+d  curve).  The  water  molecules  very  effectively  screen 
the  field  inside  the  capacitor,  except  for  the  region  z  >  0.68  run  where  the  water  distribution 
drops  rapidly  to  zero.  Adding  the  quadrupole  (also  from  water  only)  term  brings  the  atomic  and 
molecular  calculations  into  some  measure  of  agreement  (note  that  the  atomic  and  m+d+q  po¬ 
tentials  are  off  set  by  0.04  in  Figure  3).  The  quadrupole  contribution  to  the  potential  is  a  neg¬ 
ative  constant  in  the  region  away  from  the  wall  potential.  This  occurs  because  its  contribution 
to  the  charge  <p„.„((r,/)>  contains  a  double  derivative  in  space  coordinates.  The  differences  are 
greatest  at  the  surfaces  where  tlie  atomic  method  traces  charge  density  smoothly  whereas  the 
molecular  method,  based  on  an  expansion  about  molecular  centers,  requires  many  high  multi¬ 
poles  to  describe  the  field. 

The  variation  in  potential  is  rapid  near  the  metal  and  Figure  3  shows  tliis  variation  for  the  atomic 
and  m+d+q  approximation  scaled  by  a  factor  0.1.  After  adjusting  for  the  offset  of  0.04  the  two 
curves  become  identical  for  z  >  0.8  nm  as  expected  for  a  region  containing  no  molecules  or  ions. 
Between  0.6  and  0.75nm  the  difference  between  the  methods  are  largest  due  the  difference  in 
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source  terms.  The  dangers  of  omitting  higher  multipole  was  clearly  pointed  out  by  Wilson, 
Pororille  and  Pratt for  water  without  ions. 

Finally  we  return  to  the  asymmetry  in  the  lithium  ion  distribution.  There  are  effects  due  to 
system  size  and  ion-ion  correlation.  We  have  repeated  some  of  the  calculations  with  two  lithium 
ions  and  598  water_ST2  molecules  for  a  few  hundreds  of  picoseconds.  The  water  structures 
appear  almost  exactly  as  in  the  smaller  simulations.  The  lithium  probability  p(z)  for  the  two  ion 
system  is  spread  across  most  of  the  cell  in  a  diffuse  zone  with  a  moderate  bias  towards  the  metal 
side  of  the  cell.  This  result  is  satisfying  in  that  no  peculiar  or  pathological  features  are  revealed 
implying  that  the  smaller  calculations  can  give  usefully  insight  into  water  structure  within  about 
two  water  layers  (ca.  0.6  nm  from  the  point  where  the  wall  potential  becomes  repulsive  at  z=0.68 
run.). 

We  close  this  section  with  a  discussion  of  changing  the  length  scale  over  which  the  test  function 
averaging  is  performed.  The  lithium  system  is  chosen  because  tlie  ion  occupies  a  diffuse  region 
and  does  not  contact  adsorb  on  the  electrode.  The  initial  smearing  of  ion  charge  by  the  localized 
gaussian  average  does  not  overlap  the  metal  surface  region  so  that  changes  near  the  metal  reflect 
averaging  the  local  distributions.  Figure  4  shows  the  result  of  changing  the  size  of  the  length 
scale  when  gaussian  averaging  is  performed.  It  shows  that  as  the  averaging  length  scale  in¬ 
creases  the  atom  potential  loses  features  at  the  scale  of  a  water  molecule.  The  parameter  g  (in 
nm  units)  is  the  width  of  the  function /,(z)  defined  in  Eqn.(3.8).  There  are  two  families  of  curves, 
the  broken  lines  are  for  the  atom  (method  I)  potential  with  g  =  .03,  .1,  .3  and  1.0.  Though  hardly 
apparant  because  of  the  factor  x20  in  ordinate  scale,  the  atom  potential  witli  g  =  .03  is  close  in 
numerical  value  to  the  atomic  potential  shown  previously  by  the  dotted  curve  in  Figure  3.  The 
g  =  1 .0  is  extreme  because  it  is  close  to  half  the  box  edge  ViL  =  .93 1  nm.  The  behaviour  of  the 
g  =  1.0  monopole  shows  there  is  some  smeared  charge  beyond  tlie  left  hand  wall.  What  Figure 
4  demonstrates  is  tliat  as  the  size  of  the  region  averaged  is  increased  from  g  =  .03  nm  (.1  x  water 
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molecule  dimension)  to  g  =  1.0  nm  (  3  x  water  molecule  dimension)  the  structure  in  the  electric 
potential  near  the  surface  due  to  the  water  profile  is  lost.  The  total  potential  becomes  monotonic 
and  resembles  the  family  of  monopole  potential  curves  calculated  using  the  monopole  charge 
distribution.  Since  the  component  of  the  water  dipole  perpendicular  to  the  surface  is  not  aver¬ 
aged  to  zero  the  atom  potential  at  the  surface  remains  smaller  in  magnitude.  The  macroscopic 
dielectric  polarization  remains  and  reduces  the  value  of  the  surface  potential  by  roughly  an  order 
of  magnitude  (curves  for  g  =  .3  nm).  As  expected  the  monopole  potentials  are  not  as  sensitive 
to  the  value  of  the  width  g  as  the  atomic  potentials.  This  is  an  important  qualitative  result  be¬ 
cause  it  shows  the  transition  from  the  microscopic  scale,  where  surface  water  oscillatory  struc¬ 
ture  dominates  the  potential,  to  macroscopic  scale  behaviour  where  water  contributes  a  simple 
scaling  of  the  electric  potential. 


22 


V.  IODIDE  ION 


In  this  section  we  describe  simulations  for  one  iodide  anion  I  and  157  water  molecules  inter¬ 
acting  with  a  metal  surface.  In  practise  it  was  found  not  necessary  to  run  the  simulation  beyond 
500  ps  because  the  ion  adsorbed  on  the  metal  after  approximately  10  ps  and  remained  there. 
However  for  consistency  with  the  other  simulations  reported  in  this  paper  the  calculation  was 
allowed  to  run  for  1000  ps.  As  in  the  case  of  lithium  the  first  100  ps  were  used  to  allow  the 
system  to  equilibrate  and  were  excluded  from  any  statistical  analyses.  When  the  ions  and  water 
are  described  by  the  Heinzinger  parameter  set^^  we  have  shown  in  previous  publications  that: 
i)  in  fields  of  one  GV/m  all  the  halide  ions  except  fluoride  contact  adsorbed  on  non-metal  (sur¬ 
face  with  no  electrostatic  image  interactions)  electrodes^^  and  ii)  that  similar  behaviour  occurred 

for  the  model  of  the  metal  used  here  when  there  was  an  unscreened  external  field  of  1  GV/m 
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pulling  the  ion  to  the  surface  . 

Figures  5  and  6  show  the  probability  density  distributions  and  electric  potential  for  the  iodide 
system.  Figure  5  shows  the  iodide  ion  distribution  peaked  at  ca.  0.7  nm  closer  to  the  electrode 
than  where  the  wall  potential  becomes  positive  at  0.68  nm.  Compared  to  the  previous  example 
(lithium  ion  Figure  2)  there  is  less  structure  in  the  water  and  in  the  H_ST2  and  PC_ST2  com¬ 
ponent  distributions  (not  shown).  In  fact  apart  from  a  small  peak  at  0.65  nm  the  distribution 
resembles  water  without  ions  in  zero  applied  field^^’  The  atomic  charge  distribution  displays 
a  medium  positive  peak  at  ca.  0.65  nm  (H_ST2)  and  two  negative  peaks  at  ca.  0.7  nm  (iodide), 
and  ca.  0.75  nm  (PC_ST2). 

The  Figure  6  shows  Uic  electric  potential  calculated  by  the  two  methods.  All  components  go 
off  scale  near  z  =  0.682  nm,  and  are  plotted  for  z  >  0.65  nm  on  a  reduced  scale  (xO.2  in  the  right 
hand  panel).  The  monopole  potential  m  (broken  curve)  due  to  tlie  iodide  alone,  is  zero  except 
close  to  the  wall  where  the  iodide  is  localized.  There  it  changes  swiftly  from  0  to  a  high  positive 
value  corresponding  to  an  electric  field  of  approximately  -5  GV/m.  The  absense  of  an 
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electric  field  due  to  iodide  ions  means  there  is  little  or  no  reaction  from  the  water  molecules  and 
they  are  less  ordered  between  -0.6  nm  and  0.4  nm,  and  what  ordering  there  is  comes  from  the 
presence  of  the  surfaces.  In  Figure  6  the  m+d  curve  for  combined  monopole  and  dipole  poten¬ 
tials  shows  a  weaker  variation  compared  to  lithium  ion.  This  means  that  before  z  =  0.4  nm  the 
charge  on  the  electrode  is  completely  shielded.  The  curves  of  the  atomic  and  m+d+q  potentials 
overlap  well  for  z  <  0.4  nm  due  primarily  to  the  water  quadrupole  contribution. 
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VI.  NEUTRAL  SOLUTIONS.  POTENTIAL  OF  ZERO  CHARGE 

In  tliis  section  we  compute  the  properties  of  a  neutral  solution  of  Lil  comprised  of  one  Li  ion 
and  one  I  ion  and  156  water  molecules  in  the  simulation  cell,  and  compare  it  to  the  case  of  LiF 
with  the  same  number  of  water  molecules.  There  is  no  net  attractive  electric  field  on  either  ion 
because  the  total  electroststic  image  charge  on  the  metal  is  zero.  The  two  systems  are  quite 
different  because  iodide  contact  adsorbs  on  the  metal  and  creates  a  surface  localized  electric 
dipole  and  field.  Because  the  aqueous  phase  is  neutral  this  simulation  models  a  system  at  the 
potential  of  zero  charge.  This  potential  does  not  have  to  be  zero,  it  will  be  very  small  if  the 
charge  distribution  in  the  aqueous  subphase  is  everywhere  almost  zero.  The  potential  of  zero 
charge  is  considered  the  natural  reference  point  Ifom  which  to  measure  electrode  potentials. 
The  measurement  of  change  of  surface  tension  of  liquid  metals  with  electrode  potential  is  the 
only  reliable  direct  method^^  of  determining  the  PZC.  The  calculation  of  electric  field  and  po¬ 
tential  using  the  atom  charge  distribution  have  been  published^  so  in  this  section  we  focus  on 
comparing  the  atom  calculations  with  the  molecular  method  to  gain  further  insight.  Figure  7 
displays  the  most  important  density  profiles  for  the  Lil  system,  and  Figure  8  shows  the  various 
components  of  the  electric  potential. 

In  Figure  8  we  note  that  the  monopole  potential  m  drops  steeply  as  it  passes  through  the  diffuse 
layer  of  lithium  ion  (much  as  it  did  in  the  case  of  one  lithium  ion)  because  the  iodide  adsorbed 
on  the  metal  acts  like  the  second  plate  in  a  capacitor.  The  m+d  curve  shows  that  the  water 
dipoles  orient  in  the  ion  field  and  shield  tlie  ions  from  each  other  accept  close  to  the  metal. 
As  in  the  previous  two  examples  tire  water  quadrupole  shifts  die  electric  potential  to  lower  values 
by  a  constant  amount.  Note  that  the  shift  in  potential  from  vacuum  to  metal,  which  measures 
the  potential  of  zero  charge,  depends  on  Uae  net  dipole  moment  of  the  system^^.  The  atomic 
and  molecular  methods  give  similar  results.  Closer  results  requires  the  inclusion  of  octopole  and 
higher  moments  which  contribute  most  near  the  surface. 
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We  have  also  performed  simulations  on  LiF  under  the  same  conditions  that  mimic  the  potential 
of  zero  charge.  Figure  9  shows  the  probability  distributions.  In  neutral  LiF  solution  in  contrast 
to  Lil  solution  neither  ion  adsorbs  on  the  electrode  and  both  are  comingled  into  a  diffuse  layer 
that  is  predominantly  neutral  across  the  system.  The  water  distribution  looks  like  water  between 
uncharged  plates.  Figure  10  displays  the  electric  potentials.  As  expected  the  contact  potential 
was  calculated  to  close  to  zero  since  the  net  dipole  moment  is  very  small.  The  monopole  po¬ 
tential  m  shows  a  minimum  related  to  the  apparant  bimodal  distribution  of  the  lithium  ion.  The 
water  dipole  orients  to  reduce  the  local  electric  field  to  a  small  value,  with  the  result  the  curve 
m-Hd  shows  weak  changes  across  the  film.  Adding  the  quadrupole  in  m-i-d+q  shifts  the  whole 
potential  downwards  bringing  it  into  close  correspondence  with  the  result  of  the  atomic  method 
of  calculation.  For  z  <  0.4  nm  the  atomic  and  m-rd+q  potentials  are  similar  in  shape  and  value. 
Higher  multipole  contribution  to  the  molecular  method  are  needed  to  bring  it  into  closer  corre¬ 
spondence  to  the  atom  method  near  the  metal  surface. 
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VII.  LITHIUM  CATION  'COADSORBED'  WITH  IODIDE  IONS. 

In  this  simulation  we  explore  another  important  aspect  of  the  adsorption  of  ions  on  metal 
electrodes,  namely  the  ability  of  strong  contact  adsorbers  like  iodide  ions  I  to  adsoib  on  posi¬ 
tively  charged  electrodes  in  sufficent  excess  to  change  the  sign  of  the  charge  at  the  interface  as 
observed  by  an  ion  located  in  the  diffuse  layer.  In  the  case  under  consideration  cations  will  be 
attracted  out  of  the  diffuse  layer  to  compensate  the  excess  negative  charge  at  the  interface.  The 
interface  can  develope  a  layered  stmcture  of  alternating  charges.  First  there  is  positively  charged 
metal,  then  a  narrow  layer  of  contact  adsorbed  negative  ions,  and  then  a  layer  of  compensating 
cations.  In  real  systems  this  compensating  layer  may  be  part  of  a  diffuse  layer  or  it  could  be  a 
separate  structure.  The  coadsoiption  of  the  cations  is  required  to  maintain  overall  charge  neu¬ 
trality.  It  is  in  this  sense  that  we  refer  to  the  cations  as  coadsoibed.  Phenomena  of  this  type 
occur  frequently  in  electrochemistry.  In  a  previous  paper  we  very  briefly  described  results  for 
this  system  in  an  external  applied  field  corresponding  to  uncompensated  charge  on  the  immersed 
electrode.  In  the  present  calculation  the  charge  on  the  electrode  due  to  electrostaic  images  is 
+lel.  The  presence  of  two  adsorbed  iodides  creates  an  anisotropic  surface  electric  field  and  some 
simulations  were  mn  for  3  ns  to  improve  the  cation  statistics.  Even  though  only  small  changes 
in  the  distribution  were  observed  due  to  configurations  accumulated  at  later  times,  we  regard  the 
calculated  distributions  for  the  cation  as  being  approximate.  A  prelimary  account  of  the  proba¬ 
bility  distributions  of  this  simulation  has  been  published^^.  This  is  the  first  report  of  the  po¬ 
tentials  across  this  cell. 

Figure  1 1  shows  tlie  probability  density  profiles  for  selected  species  in  solution.  There  are  some 
similarities  in  the  results  for  water  and  tlie  lithium  ion  for  this  system  and  the  one  discussed  in 
the  last  section.  In  particular  the  water  looks  quite  unstructured  compared  to  hydrated  ions,  and 
the  Li’'  ion  is  bimodal.  The  L  ion  distribution  is  sharper  because  the  electrode  is  positively 

charged. 
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The  initial  configuration  was  a  random  arrangement  of  ions  and  water  on  a  cubic  lattice.  Both 
the  iodides  adsorbed  within  50  ps.  Note  that  the  iodide  distribution  is  localized  at  the  surface 
as  in  the  case  of  one  iodide,  and  the  lithium  is  more  diffuse  than  in  the  case  of  a  single  lithium 
in  the  field  of  its  own  electrostatic  image  charge.  Except  for  slightly  accentuated  broad  struc¬ 
tures  in  the  range  0.4  to  0.6  nm,  the  water  distribution  resembles  the  distribution  already  calcu¬ 
lated  for  zero  applied  field^^.  In  particular  we  note  that  there  are  no  new  peaks  near  0.65  nm 
indicative  of  localized  surface  water.  The  proton  H_ST2  and  point  charge  PC_ST2  components 
have  more  structure  near  the  metal  than  in  zero  field.  This  is  to  be  expected  since  they  are 
charged  and  the  field  between  them  and  their  images  is  unscreened.  The  atomic  charge  density 
has  a  deep  minimum  at  ca.  0.7  nm  due  to  iodide  and  a  shoulder  at  ca.  0.75  nm  due  to  the 
PC_ST2. 

Figure  12  shows  the  electric  potential  calculated  using  the  two  methods.  Starting  from  the  left 
the  monopole  potential  first  rises  as  it  passes  through  the  lithium  layer,  and  then  drops  rapidly 
as  it  approaches  the  sheet  of  adsorbed  iodide,  finally  it  turns  upward  because  the  electrode  carries 
a  net  positive  image  charge.  Apart  from  the  upward  turn  near  the  metal  this  monopole  behaviour 
is  qualitatively  similar  to  that  depicted  for  the  model  in  Figure  1.  Our  assertion  that  the  tradi¬ 
tional  models  neglect  significant  contributions  from  the  water  can  be  seen  by  examination  of 
Figure  12.  There  is  significant  structure  in  the  electric  potential  calculated  by  either  method 
induced  by  orientation  and  packing  distribution  of  the  waters  at  the  surface.  The  monopole  term 
is  completely  compensated  (-0.65  <  z  <  0.4  nm)  by  the  dipole  potential  from  the  water  molecules 
as  shown  by  curve  m+d  until  next  to  the  metal  surface  where  water  is  more  strongly  oriented 
and  where  the  water  supply  drops  rapidly  to  zero.  From  0.4  nm  to  0.65  nm  die  atomic  potential 
closely  follows  tliat  of  one  iodide  (see  Figure  6).  On  the  right  hand  side  of  Figure  12  we  show 
the  potentials  atomic  (atm),  monopole  m  and  m+d-t-q  scaled  by  a  factor  0.1. 
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Scale  up  is  likely  to  be  an  important  issue  since  the  size  of  the  system  may  be  too  small  to  allow 
the  diffuse  layer  lithium  ion  enough  space  to  distribute  without  interference  from  the  restraining 
wall  at  z=-0.931  nm.  However  conclusions  drawn  concerning  the  behaviour  of  ions  within  1.0 
ran  of  the  metal  surface  will  not  be  very  sensitive  to  the  position  of  the  restraining  wall.  An 
example  of  a  large  system  with  intact  double  layer  is  discussed  next. 
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VIIL  IM  NaCl  SOLUTION. 


In  this  final  example  we  describe  a  simulation  with  approximately  1600  water  molecules  and 
ions  in  a  box  with  edge  length  L  =  3.724  nm.  The  simulation  cell  contains  32  Na  ions,  28 
Cf  ions,  and  1576  water  molecules.  The  NaCl  concentration  is  about  IM.  The  electrostatic 
image  charge  on  the  electrode  surface  due  to  the  difference  in  number  of  positive  over  negative 
ions  is  -4lel  or  0.046  C  m'^.  This  essentially  the  same  charge  density  as  in  all  the  earlier  smaller 
simulations  with  non  zero  electrode  charge  density.  Computer  time  constraints  dictated  the  use 
of  a  simple  water  model.  For  this  reason  we  chose  the  SPCE  water  model  (three  charged  mass 
points^^)  over  the  ST2  water  model  (two  charged  mass  points,  two  charged  zero  mass  points, 
one  uncharged  mass  point).  The  NaCl  parameters'^  used  are  those  appropriate  for  SPCE  water. 
We  also  experimented  with  simulations  nm  at  elevated  temperatures  to  increase  the  diffusion  rate 
of  the  ions  in  order  to  get  better  statistics.  In  this  calculation  the  temperature  was  30  °C.  The 
simulation  was  run  for  840  ps  with  the  first  100  ps  used  to  equilibrate  the  system.  On  a  dedi¬ 
cated  IBM  RS6000  model  550  workstation  the  calculation  takes  approximately  10  to  12  weeks. 

Figure  13  shows  the  probability  density  profiles  for  the  water  proton,  water  mass  center,  Na^ 
ions  and  Cl  ions.  The  anion  distribution  (broken  line)  has  been  smoothed  to  permit  it  to  be 
distinguished  from  the  cation.  Also  shown  rising  monotonically  from  the  left  are  the  integrals 
of  the  ion  densities.  The  near  coincidence  of  the  integrals  for  z  <  0.7  nm  shows  that  the 
electrolyte  is  approximately  charge  neutral.  For  z  >  0.7  nm  the  integrated  densities  systemat¬ 
ically  diverge  as  expected  for  a  transition  from  the  locally  neutral  'bulk'  electrolyte  into  the 
'diffuse'  part  of  the  electric  double  layer  where  screening  occurs.  At  z  =  1.0  nm  the  divergence 
in  the  integrated  densities  equals  the  largest  previous  difference  in  the  two  curves  implying  that 
the  region  z  >  1.0  nm  corresponds  to  solution  shielding  the  charge  on  the  metal.  The  Na”"  ion 
distribution  shows  well  defined  structure  in  the  form  of  a  broad  peak  at  ca.  1.1  nm,  and  a  sharp 
peak  at  1.4  nm.  The  water  and  proton  distributions  appear  flat  for  Izl  <  0.8  nm.  On  the  metal 


30 


side  the  water  probability  distribution  has  peaks  at  0.9  nm,  1.2  nm,  and  a  strong  asymetric  peak 
at  1.6  nm.  This  latter  feature  appears  to  be  composite  being  the  superposition  of  a  broad  feature 
at  1.5  nm  and  a  narrow  peak  at  1.6  nm.  The  peaks  in  the  proton  distribution  are  at  0.9  nm,  1.2 
nm,  1.55  nm,  and  1.7  nm.  The  last  peak  at  ca.  1.7  nm  comes  from  protons  in  water  OH  bonds 

pointing  at  the  metal. 

Figure  14  shows  the  atomic  charge  density  for  a  bin  function  with  width  0.004655  nm  (  or 
L/800).  We  note  that  the  charge  density  appears  flat  for  Izl  <  0.8  nm.  The  contribution  from 
the  ionic  charge  for  z  >  0.8  nm  is  not  evident  because  the  charge  on  the  waters  dominates.  The 
electric  field  was  obtained  by  integration  from  to  position  z.  The  field  is  flat  with  small 
variations  around  zero  in  the  region  Izl  <0.8  nm.  Near  the  metal  the  electric  field  undergoes  a 
series  of  rapid  oscillations  as  it  vears  upwards.  These  oscillations  are  due  to  the  water  stracture 
at  the  interface.  The  general  rise  is  due  to  the  excess  Na^  charge  in  the  double  layer. 

Figure  15  shows  the  potential  calculated  using  the  atom  method,  and  the  components  of  the 
potential  calculated  by  the  molecule  method.  The  contact  potential  is  about  -2.0  V,  and  as  in 
the  smaller  simulations  the  potential  in  the  'bulk-like'  zone  comes  from  the  water  quadrupole. 
The  broken  curve  labeled  m  is  that  from  the  ionic  charge.  If  the  system  were  truely  neutral  then 
curve  m  would  be  flat  and  zero  all  the  way  to  the  beginning  of  the  diffuse  layer.  Clearly  the 
small  charge  imbalances  seen  in  Figure  13  do  affect  shape  of  the  potential  curve  and  may  be 
an  exacting  criterion  with  which  to  judge  these  larger  simulations.  The  transition  to  monotomc 
decrease  starting  near  z  =  0.7  nm  is  another  indicator  of  where  the  diffuse  layer  begins  in  this 
simulation.  Note  that  m+d  shows  that  the  dipole  potential  completely  compensates  the 
monopole.  Including  the  quadrupole  to  give  m+d+q  shifts  the  core  region  downwards  by  3  V 
and  brings  the  molecular  calculation  of  potential  into  correspondence  with  the  atom  metliod  of 
calculation.  The  molecular  method  has  larger  extrema  near  tire  surface  compared  to  the  atom 
method.  The  reason  for  this  is  tlie  very  small  bin  size  and  the  consequential  need  to  include 


31 


many  high  order  multipoles  in  the  molecule  method.  However  as  required  the  contact  potential 
is  the  same  in  each  case.  Note  also  that  adding  the  quadrupole  term  does  not  change  the  contact 
potential. 

In  summary  this  larger  simulation  has  verified  that  water  structure  near  walls  is  not  an  artifact 
of  small  size.  In  IM  NaCI  the  double  layer  is  about  1  ran  thick  and  encompasses  about  three 
layers  of  water.  These  water  layers  can  significantly  affect  the  distribution  of  ions  near  the  metal 
creating  features  in  the  probability  distributions  that  are  not  describable  in  the  Gouy-Chapman- 
Stem  model. 
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IX.  CONCLUSIONS 


In  this  paper  we  have  shown  how  the  results  of  molecular  dynamics  simulations  performed  at 
constant  N,V,T  can  be  used  to  calculate  the  electrostatic  field  and  potential  across  an  electric 
double  layer.  Apart  from  usual  molecular  dynamics  assumptions  we  used  a  model  for  an  im- 
mereed  electrode  in  which  electrode  charge  was  equal  to  the  image  charge  on  the  metal.  All 
electrostatics  were  summed  without  trunkating  the  long  range  coulomb  interaction  using  the  fast 
multipole  method.  The  charge  distribution  in  the  system  was  analysed  by  two  methods.  One 
was  based  on  the  charges  in  atoms  in  vacuum,  the  other  was  based  on  molecules  and  required 
expanding  the  positions  of  charge  inside  a  molecule  about  an  origin  on  each  molecule.  The  later 
method  provided  a  natural  route  for  decomposing  fields  into  monopole  and  higher  order  multi¬ 
pole  components.  We  also  showed,  by  increasing  the  length  scale  over  which  local  average  is 
taken,  how  microscopic  structure  due  to  the  water  probability  distribution  is  systematically 
washed  away  leaving  only  the  macroscopic  effect  of  the  water  dielectric  polarization  on  the 
electric  surface  potential.  We  also  showed  using  the  results  of  a  calculation  almost  an  order  of 
magnitude  larger  in  number  of  water  molecules,  that  the  main  features  evident  in  smaller  cal¬ 
culations  of  the  probability  profiles  for  water  near  charged  metal  surfaces  are  not  greatly  affected 
by  scale  up.  This  implies  that  small  simulations  involving  a  few  hundred  water  molecules  pro¬ 
vide  a  useful  starting  point  for  studying  double  layer  properties. 

In  a  series  of  examples  starting  with  single  ions  in  water,  and  tlien  to  increasingly  more  complex 
cases  we  explored  fields  close  to  the  surface  and  across  double  layers  in  some  electrochemically 
interesting  examples.  The  case  of  lithium  was  chosen  to  explore  some  of  the  consequences  of 
averaging  over  length  scales  from  smaller  to  larger  than  a  water  molecule.  In  die  final  example 
of  IM  NaCl  we  demonstrated  the  coexistence  of  bulk  and  double  layer  regions  in  equilibrium 
in  the  same  simulation.  Since  the  system  contained  sixty  ions  there  is  a  ample  correlation  be¬ 


tween  ions  and  amongst  ion  and  water  molecules. 
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Figure  1.  Schematic  diagram  of  the  traditional  view  of  the  electric  double  layer  showing  a  flat 
negatively  charged  electrode  with  contact  adsorbed  anion  and  fully  hydrated  cations  at  the  outer 
Helmholtz  plane  (OHP).  Electric  potential  shown  schematically  at  the  top. 
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Figure  2.  Density  profiles  for  lithium  cation  Li"^  ion,  157  ST2  waters  and  the  total  atomic  charge 
density  p  near  an  immersed  electrode.  Image  charge  on  metal  -1  is  screened  by  the  Li"^  ion. 
Metal  electrode  on  right  hand  side,  restraining  wall  on  the  left.  Image  plane  at  z  =  0.931  nm. 
Wall  potentials  go  through  zero  at  171=0.682  nm.  Simulation  from  100  ps  to  2500  ps. 
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Figure  3.  Potential  drop  across  the  system  Li+  ion  and  157  ST2  waters  in  the  immersed  electrode 
model.  Calculations  using  atom  and  molecule  methods.  Bin  size  0.004655  nm.  Molecule 
method;  m  monopole  only,  m+d  monopole  and  dipole,  m+d+q  monopole  dipole  and  quadrupole. 
The  right  side  panel  shows  the  potentials  on  larger  scale  for  z  >  0.4  nm,  with  the  monopole 
scaled  by  0.1.  Image  charge  on  metal  -1  is  screened  by  the  Li+  ion.  Metal  electrode  on  right 
hand  side,  non  metallic  restraining  wall  on  the  left.  Image  plane  at  z  =  0.931  nm.  Wall  po¬ 
tentials  go  through  zero  at  lzl=0.682  nm. 


.4  .6  .8 


4i 


Potential 


Distance  z  /  nm 


Figure  4.  Potential  drop  across  the  system  Li+  ion  and  157  ST2  waters  in  the  immersed  electrode 
model.  Calculations  using  atom  and  molecule  methods  with  gaussian  widths  g  =  .03,  .1,  .3,  1.0 
nm.  Curves  for  monopole  potentials  (broken  lines)  and  the  total  potential  by  the  atom  (dotted 
lines)  method.  Oscillations  in  the  atomic  potential  clearly  visible  in  Figure  3  for  z  >  0.2  nm  are 
washed  out  on  the  scale  of  a  water  molecule  g  =  .3. 
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Figure  6.  Potential  drop  across  the  system  r  ion  and  157  ST2  waters  in  the  immersed  electrode 
model.  Calculations  using  atom  and  molecule  methods.  Bin  size  0.004655  nm.  Molecule 
method:  m  monopole  only,  m+d  monopole  and  dipole,  m+d+q  monopole,  dipole  and  quadrupole. 
Panel  on  the  right  hand  side  shows  the  curve  for  z  >  0.4  nm  on  a  reduced  scale.  Image  charge 
on  metal  +1  is  screened  by  the  F  ion.  Metal  electrode  on  right  hand  side,  restraining  wall  on 
the  left.  Image  plane  at  z  =  0.931  nm.  Wall  potentials  go  through  zero  at  lzl=0.682  nm. 
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Figure  7.  Density  profiles  a  neutral  solution  consisting  of  one  lithium  cation  Li-",  iodide  anion 
F,  and  157  ST2  waters  and  the  total  atomic  charge  density  near  an  immersed  electrode.  Metal 
electrode  on  right  hand  side,  restraining  wall  on  the  left.  Image  plane  at  z  =  0.931  nm.  Wall 
potentials  go  tlirough  zero  at  lzl=0.682  nm. 
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Figure  8.  Electric  potential  drop  across  a  the  neutral  solution  consisting  of  one  lithium  ion 
Li'^,  one  iodide  ion  F,  and  156  ST2  waters  next  to  an  immersed  electrode.  Potential  due  to  atom 
method  shifted  by  -0.04.  Potential  from:  ions  only  m,  ions  and  water  dipoles  m+d,  ions  and 
water  dipole  and  quadrupole  m+d+q.  Note  that  the  potential  m+d+q  is  almost  the  same  as  cal¬ 
culated  from  atomic  charges.  There  is  no  net  image  charge  on  the  metal.  Metal  electrode  on 
right  hand  side,  restraining  wall  on  the  left.  Image  plane  at  z  =  0.931  nm.  Wall  potentials  go 
through  zero  at  lzl=0.682  nm. 
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Figure  10,  Electric  potential  drop  across  a  neutral  solution  consisting  of  one  lithium  ion  Li"^, 
one  fluoride  ion  F",  and  156  ST2  waters  next  to  an  immersed  electrode.  Potential  due  to  atom 
method  is  shifted  by  -0.04.  Potential  from:  ions  only  m,  ions  and  water  dipoles  m+d,  ions  and 
water  dipole  and  quadrupolc  m-i-d+q.  Note  that  the  potential  m-i-d+q  is  almost  the  same  as  cal¬ 
culated  from  atomic  charges.  Metal  electrode  on  right  hand  side,  restraining  wall  on  the  left. 
Image  plane  at  z  =  0.931  nm.  Wall  potential  goes  tlirough  zero  at  lzl=0.682  nm. 


-0.8  -0.4  0.0  0.4  0.8 

Distance  z  /  nm 


48 


Number  Density  px(z)  /  nm 


Potential  V  /  5.22 


-.8  -.4  0  .4  .6  .8  .4  .6  .8 

Distance  z  /  nm 


Figure  12.  Electric  potential  drop  across  a  charged  solution  consisting  of  one  lithium  ion  Li"^, 
two  iodide  ions  T,  and  155  ST2  waters  next  to  an  immersed  electrode.  Potential  due  to  atom 
method  is  shifted  by  -0.04.  Potential  from:  ions  only  m.  ions  and  water  dipoles  m-td,  ions  and 
water  dipole  and  quadrupole  m-td+q.  Note  that  the  potential  m-i-d+q  is  almost  the  same  as  cal¬ 
culated  from  atomic  charges.  Met^  electrode  on  right  hand  side,  restraining  wall  on  the  left. 
Image  plane  at  z  =  0.931  nm.  Wall  potential  goes  through  zero  at  lzl=0.682  nm.  Panel  on  right 
shows  potentials  for  z  >  0.4  nm  on  a  scale  x5  larger  than  the  left  side. 
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Figure  14.  Vacuum  charge  density,  electric  field  and  potential  for  IM  NaClaq  at  30  °C  and  -4e 
image  charges.  Electrolyte  composition:  32  Na'*'  ions,  28  Cl"  ions,  and  1576  SPCE  water  mol¬ 
ecules  near  an  immersed  electrode.  Metal  electrode  on  right  hand  side  at  z  =  1.862nm  (position 
of  image  plane),  restraining  wall  origin  on  the  left  at  z  =  -1,862  nm,  wall  potentials  on  both  sides 
of  the  simulation  cell  go  through  zero  at  lz>  1.615  nm.  Simulation  duration  100  ps  to  840  ps. 
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Fif’ure  15.  Electric  potential  drop  across  a  charged  solution  consisting  of  for  IM  NaCl  solution 
at  Iso  °C  and  -4lel  image  charge  on  tlie  metal.  Potential  calculated  hy  the  atom  method  is  shitted 
by  -0.8.  Molecule  method  potentials  from:  ions  only  m,  ions  and  water  dipoles  m-i-d,  ions  and 
water  dipole  and  quadrupole  m-td-t-q.  Note  that  the  potential  m-i-d-tq  is  almost  the  same  as  cal¬ 
culated  by  tlie  atom  method.  The  net  image  charge  on  tlie  metal  -4lcl  is  screened  hy  a  total  ot 
sixty  ions.  Electrolyte  composition:  32  Na"^  ions,  28  CP  ions,  ^d  1576  SPCE  water  molecules 
near  an  immersed  electrode.  Metal  electrode  on  right  hand  side  at  z  =  1.862nm  (position  or 
image  plane),  restraining  wall  origin  on  the  left  at  z  =  -1.862  nm,  wall  potentials  on  bo^  sides 
of  the  simulation  cell  go  through  zero  at  lzl=  1.615  nm.  Simulation  duration  100  ps  to  840  ps. 


